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An ensemble of pulse-coupled phase-oscillators is thoroughly analysed in the presence of a mean-field cou¬ 
pling and a dispersion of their natural frequencies. In spite of the analogies with the Kuramoto setup, a much 
richer scenario is observed. The “synchronized” phase, which emerges upon increasing the coupling strength, 
is characterized by highly-irregular fluctuations: a time-series analysis reveals that the dynamics of the order 
parameter is indeed high-dimensional. The complex dynamics appears to be the result of the non-perturbative 
action of a suitably shaped phase-response curve. Such mechanism differs from the often invoked balance be¬ 
tween excitation and inhibition and might provide an alternative basis to account for the self-sustained brain 
activity in the resting state. The potential interest of this dynamical regime is further strengthened by its (mi¬ 
croscopic) linear stability, which makes it quite suited for computational tasks. The overall study has been 
performed by combining analytical and numerical studies, starting from the linear stability analysis of the asyn¬ 
chronous regime, to include the Fourier analysis of the Kuramoto order parameter, the computation of various 
types of Lyapunov exponents, and a microscopic study of the inter-spike intervals. 
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I. INTRODUCTION 

Most of the challenging questions that arise in the attempt 
of improving our understanding of the natural (and artificial) 
world deal with multi-component systems, whose overall dy¬ 
namics is the result of many nonlinear interactions. The diffi¬ 
culty of the task is often mitigated by the assumption of being 
before universal phenomena, which do not crucially depend 
on the details of the underlying models. It is therefore cus¬ 
tomary to deal with relatively simple setups in the hope that 
relevant details are not missed. 

The mammalian brain is the most prominent example 
where this approach is absolutely necessary, if we wish to 
make some substantial progress. There, even after disre¬ 
garding several ingredients, such as the multiple degrees of 
freedom involved in the dynamics of realistic neurons (as in 
multicompartmental models [1]), the diversity among the sin¬ 
gle units, the topology and the plasticity of the connections, 
the range of possible dynamical phenomena is still very rich 
and not yet entirely understood. Self-consistent partial syn¬ 
chronization is a simple but enlightening example. The phe¬ 
nomenon, discovered by van Vreeswijk in an ensemble of 
leaky integrate-and-fire neurons (LIF) [2], was believed for 
a long time to be a non-perturbative effect. Only, recently it 
has been however clarified [3] that it is equivalent to the ro¬ 
tating waves observed in the weak-coupling limit [4], and can 
indeed be observed and characterized in Kuramoto-Daido os¬ 
cillators [5], as well. 

In general, the problem of characterizing the collective dy¬ 
namics of an ensemble of oscillators is deeply connected to 
the question of how different levels of descriptions are linked 
to one another. In computational neuroscience, it is customary 
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to consider subpopulations of neurons, under the assumption 
that the firing rate is the single relevant variable, as in the sem¬ 
inal paper by Wilson and Cowan [6] and in several other pub¬ 
lications (see e.g. [1, 7, 8]). It is not, however, clear whether 
such models can be derived starting from more microscopic 
setups based on single spiking neurons. Some recent stud¬ 
ies have shown that a low-dimensional collective dynamics 
may emerge in networks of theta (or, equivalently, quadratic- 
integrate-and-fire, QIF) neurons [9, 10]. More than that, a 
reformulation of pulse-coupled oscillators in terms of firing- 
rate models has been accomplished in [11, 12]. The validity 
of these results is due to the existence of relationships such 
as the Ott-Antonsen Ansatz [13] and the Watanabe-Strogatz 
theorem [14], which allow expressing the collective behav¬ 
ior in terms of a few variables, the others being essentially 
slaved. Such theoretical pillars are however based on strong 
simplifying assumptions on the nature of the inter-oscillator 
coupling [13, 14]. 

To what extent is the compression of degrees of freedom 
effective in more general setups? The background activity 
of the brain in the resting state, when no specific task is 
performed [15-17] testifies to a collective irregular dynam¬ 
ics. Moreover, the ongoing discussion about rate- versus 
temporal-coding [18, 19], suggests that the firing rate may 
not be sufficient to ensure the necessary computational capa¬ 
bility of the mammalian brain. 

Altogether, one should thus expect an irregular collective 
behavior. It is often conjectured that the self-sustained activ¬ 
ity is the result of a balance between activation and inhibi¬ 
tion [20, 21]. Mathematically, this means that the effect of the 
coupling is zero on average so that it is essentially controlled 
by stochastic/chaotic fluctuations. It is not, however, clear 
how such a balance can be durably ensured in self-organized 
networks of bring oscillators. A conceptually different possi¬ 
bility to account for a macroscopic irregularity is offered by 
the nonlinear character of the Fiouville-type equation (which, 
strictly speaking, applies to an ensemble of infinitely many 



2 


oscillators). This functional equation operates in an infinite¬ 
dimensional phase space and can, thereby, generate a dynam¬ 
ics of arbitrary complexity. This has been indeed observed in 
an abstract model of coupled maps [22] and in globally cou¬ 
pled Stuart-Landau oscillators [23]. In both cases, the single 
dynamical units are intrinsically more complex than phase os¬ 
cillators; the logistic maps are chaotic by themselves, Stuart- 
Landau oscillators can behave chaotically under the action of 
a periodic modulation. In the case of phase oscillators, there 
is only a preliminary evidence in an ensemble of LIF neurons 
with delayed interactions [24]. 

In this paper, we present a model whose overall activity is 
intrinsically highly-dimensional. As this dynamical phase is 
rather robust against variations of several parameters, it may 
provide an alternative mechanism for the self-sustainment of 
the resting brain activity. More precisely, we study an ensem¬ 
ble of pulse-coupled phase oscillators, whose phase-response 
curve (PRC) is derived by smoothing the PRC of LIF neurons. 
Other than that, our setup is the same as in the standard Ku- 
ramoto model [25, 26]: the single oscillators are characterized 
by a distribution of bare frequencies, while the coupling is ho¬ 
mogeneously all-to-all. As in the Kuramoto model, a synchro¬ 
nization transition is observed upon increasing the coupling 
strength, but the analogies end here, since above criticality, the 
order parameter, rather than being constant, exhibits complex 
high-dimensional oscillations. Such fluctuations are present 
in the “activity” field as well, a variable akin to the electric 
potential recorded while measuring EEGs. 

As briefly discussed in section III, in the weak-coupling 
limit (and for a small dispersion of the frequencies), our sys¬ 
tem reduces to a Kuramoto-Daido model, with the coupling 
function being composed of several Eourier harmonics. Re¬ 
cent studies of such a type of models have revealed quite a 
rich phenomenology (see, e.g. [27, 28]). This is not, however, 
sufficient to account for the qualitative differences reported 
in this paper: direct simulations show that the scenario here¬ 
after discussed disappears when the dispersion of frequencies 
is decreased. 

More specifically, the overall dynamics is characterized by 
a spectrum of negative Lyapunov exponents. This inconsis¬ 
tency is nothing but a manifestation of stable chaos [29], an 
irregular dynamics of cellular-automaton type, which is self- 
sustained because of the high (infinite) dimensionality of the 
phase space (in other words, it dies out in finite ensembles). In 
neural systems, stable chaos was first found in a diluted net¬ 
work of LIE units [30], and later discussed in more disordered 
setups [31, 32]. At variance with deterministic chaos, accom¬ 
panied by an exponential separation of orbits and thereby a 
loss of memory, stable chaos is identified by a “microscopi¬ 
cally” stable dynamics, which is definitely more appropriate 
for the performance of computational tasks. The potentiality 
of stable chaos for information processing has been prelimi¬ 
narily explored in [33, 34]. The onset of a macroscopic irreg¬ 
ular dynamics, as discussed in this paper, makes this perspec¬ 
tive even more intriguing, for the richness of the collective 
behavior. 

In section II we introduce the setup and justify its choice. 
In the following section III we reconstruct the asynchronous 


state and investigate its stability properties. Differences and 
analogies with the standard Kuramoto model are emphasized. 
In particular, we And that the asynchronous regime loses its 
stability when a complex eigenvalue is born out of the line 
containing the continuous spectrum. In the last part of the 
section, the proper order parameters for the characterization 
of the transition are introduced: they are the Kuramoto order 
parameter, whose definition requires passing to more appro¬ 
priate phase variables, and the activity field. Sections IV and 
V are devoted to a careful numerical analysis of the synchro¬ 
nized phase at the collective and microscopic level, respec¬ 
tively. Due to the difficulty of dealing with finite-size cor¬ 
rections, we study the resulting behavior sufficiently far from 
the transition. In section IV we first illustrate the phase dia¬ 
gram and the initial part of the Lyapunov spectrum. We then 
show the power spectrum of the order parameter and carry on 
a time-series analysis to determine the fractal dimension. In 
section V we focus our interest on the behavior of the sin¬ 
gle neurons, computing the effective frequency and the con¬ 
ditional Lyapunov exponents: they are all negative, indicating 
that we are in the presence of generalized synchronization. 
The presence of phase slips is also unveiled. Einally in the 
last section we summarize the main results and discuss the 
several perspectives that are opened by the scenario discussed 
in the paper. 

II. THE MODEL 

The starting point of this paper is the model of delayed LIE 
neurons studied in [24]. Here, the model is modified to make 
it simpler, more generic and more amenable to both numerical 
and analytical studies. 

In this perspective, we consider an ensemble of pulse- 
coupled phase oscillators, in the presence of (5-like pulse and 
charactereized by a suitable PRC F ((/>), 

= Wi - ^r((/)j)^(5(f-fj), (1) 
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where (pi G [0,1] is the local phase, uJi the bare oscillator 
frequency (i.e. in the absence of coupling), g the coupling 
strength and Ai is the system size. Whenever any oscillator 
reaches the threshold (pi = 1, a (5-spike is sent and received 
by all neurons. The above formulation is quite general for two 
reasons; (i) any model where the velocity held pi is phase de¬ 
pendent (in the absence of coupling), can be always rephrased 
as Eq. (1) upon suitably changing variables [35]; (ii) finite- 
width pulses can be mapped onto (5-like ones, upon suitably 
adjusting the shape of the PRC [3] (at least in the weak cou¬ 
pling limit). 

An important ingredient of the model studied in [24] is the 
presence of a delay between spike emission and reception. In 
the weak-coupling limit, when the dynamics is nearly homo¬ 
geneous, one can simulate the presence of a delay as a suitable 
phase shift of the PRC and this is what has been assumed here. 
The phase shift should be different for the different oscillators. 
However, here, for the sake of simplicity we assume the same 
PRC for all the oscillators. 
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In the LIF model, r{(j)i) = aexp{b<j)i) where (pi is assumed 
to be taken modulus 1. As a result of a phase shift, the dis¬ 
continuity originally present when passing from 1 to 0 , moves 
inside the unit interval. For the sake of generality and sim¬ 
plicity, we prefer to remove the discontinuity, considering a 
piece-wise linear PRC, such as 

f Bqi + bip if 0 <(j) < pi 
r((;i) = i Bo 2 - b2p if pi < p < pr ( 2 ) 

[-003 + bip if pr <p <l , 

where the various parameters are chosen so as to ensure con¬ 
tinuity in pi, pr and equality between ^ = 0 and 1. Con¬ 
sidering that the amplitude of the PRC is controlled by the 
coupling constant g, there are three truly independent param¬ 
eters: one controlling the vertical shift of the PRC, and two 
which identify the junction points. As for the first param¬ 
eter, it basically controls whether the coupling has an aver¬ 
age excitatory or inhibitory effect, thereby inducing a speed¬ 
ing up or slowing down of the spiking activity. Since we 
are not interested in such effects, but rather in the mutual at¬ 
traction or repulsion among the oscillators, we have decided 
to assume that the PRC has zero average. The two remain¬ 
ing parameters are identified by the phase shift s (defined as 
the distance of the midpoint of the central region from 1 - 
see the Fig. 1) and the width 6 of the central interval. Alto¬ 
gether, 62 = bi/S, Bqi = bi{s - 1 / 2 ), Bo 2 = bi{l - s)/5, 
Bo 3 = bi{s — 3/2), while pi = (1 — s + 5/2 — 6s)/{S -f 1 ) 
and pr = (1 — s + 3(5/2 — 6s)/{S -I- 1). The parameter bi 
has been set equal to 1.5 (in principle, it can be absorbed in 
the definition of g), while the two other parameters have been 
set s = 0.14, (5 = 0.1 in all of the following simulations. The 
resulting shape of the PRC is presented in Fig. 1. 



FIG. 1. Phase response curve T{pi) according to Eq. (2) for the stan¬ 
dard parameter values fci = 1.5, s = 0.14, and <5 = 0.1 (solid line); 
the short- and long-dashed curves correspond to the effective PRC F 
obtained for 10 = uimin ~ 0.8 and ui = Umax = 2, respectively 
(with g = 0.8). 

Finally, we have chosen to work with a uniform distribution 
of frequencies centered in (D = 1.4. The simulations reported 


in this paper refer to a half width A = {oJmax — i^min)/‘2 = 
0 . 6 , but similar results have been obtained for different val¬ 
ues of A. We evolve the model Eq. (1,2) as an event driven 
process. Between two consecutive (5-spikes, the phase of each 
oscillator increases linearly according to its individual bare 
frequency uJi. When one of the oscillators reaches the firing 
threshold pi = 1 , its phase is reset to zero and all phases 
are adjusted to account for the received spike. The effect of 
the coupling might bring a second oscillator beyond the firing 
threshold. In such a case, also that oscillator is reset to zero 
plus an offset due to the spike received from the first oscilla¬ 
tor. We continue this evolution, without advancing the time, 
until no further spikes are triggered. In practice, “avalanches” 
may occur: we have controlled that they do not contribute sig¬ 
nificantly to the global behavior as their size does not increase 
upon increasing the number of neurons. Notice also that in the 
original model [24], avalanches do not exist. 

III. THEORY 

Some insight can be gained by considering the thermody¬ 
namic limit, as this allows determining analytically the prop¬ 
erties of the stationary asynchronous regime. 

First of all it is convenient to define the activity field E{t) 
as the number of spikes emitted per unit time, so that Eq. (1) 
can be rewritten as (for the sake of simplicity we drop the 
subindex i), 

P = u- gV{P)E{t) . (3) 

Let us now introduce the probability density Q{p,uj,t), as 
the fraction of neurons with a bare frequency in [w, uj + duj), 
whose phase belongs to [p, p + dp) at time t. Clearly, 

J Q{p,u},t)dp = P{uj) , 

where F’(w) is the density of neurons with bare frequency ui. 
Q satisfies the continuity equation 

^ = -^[oJ-gmE{t)]Q , (4) 

while the field E satisfies the self-consistent equation 

Eit)= J Qil,uj,t)[u;-gTil)E{t)]duj , (5) 

which implies 

^ S i^Q{l,uj,t)du: 

^ + gT{l) J Q{l,uj,t)duj ■ 

The asynchronous regime corresponds to the stationary so¬ 
lution whose phase-dependence is determined by setting the 
time derivative of Q equal to zero. By properly renormalizing 
the flux, one obtains 

Qo{,p7^) 

T{l,u;)[uj - gT{p)Eo] 
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where 

is the time required by an oscillator with frequency oj to reach 
the phase if), starting from 0, in the presence of a constant 
field Eq. T{1,u!, Eq) is thereby the interspike interval, while 
T{(j)i,u!) is the inverse instantaneous effective frequency. The 
field Eq can be finally obtained from Eqs. (5,6) 


Eq = 


Pj^) 

T{l,u;) 


duj . 


( 8 ) 


The above calculation yields the structure of the asyn¬ 
chronous state, but it does not tell us whether it is stable. The 
stability can be assessed by investigating the behavior of in¬ 
finitesimal perturbations. Let us define, 

Q{(j),w,t) = Qo{(j),uj) + q{(j),uj,t) , E{t)=Eo + e{t) 


q and e satisfy the following equations. 


dq 

m 



gT{^)Eo] q + ge{t) - — - 


and 


^ j{uj - gV{l)EQ)q{l,uj,t)duj 

1-bpr(l)/Qo(l,w)dw 


These two equations can be solved by introducing a standard 
Ansatz, q{(l),u!,t) = it(^, u;)e^‘, e{t) = One obtains 

fiu = gT'EQU - [w - gTEo] u' + gV'QQZ + gTQ^z (9) 
and 

- gr(i)£^o)M(i, uj)duj 

l + 5r(l)/Qo(l,w)dw 

where the prime denotes a derivative with respect to ip and we 
have dropped the depedence on (p for the sake of simplicity. 
The solution of such equation, reported in the appendix A, 
yields the eigenvalue equation (A9). 

The spectrum of the linear operator consists of a continuous 
and a discrete component. The continuous part is confined to 
an interval along the imaginary axis and is therefore composed 
of marginally stable directions. The discrete component can 
be obtained by assuming p = gji -f ijii, separating (A9) into 
real and imaginary parts, and finally looking for the zeros in 
the complex plane. 

A numerical study reveals the presence of (at least) three 
pairs of complex-conjugate eigenvalues (see Fig. 2, where 
the imaginary part is not reported) two of which being neg¬ 
ative and one positive. Two pairs of exponents arise definitely 
above some finite g value; the third one is likely to follow 
the same scenario, but given its small real part we could not 
trace it for small coupling strengths (see the crosses in Fig. 2). 
Altogether, the asynchronous solution is marginally stable up 
until 5c ~ 0.72, when it destabilizes for the onset of a pair of 
complex conjugate eigenvalues with a positive real part. 



FIG. 2. Stability diagram of the asynchronous state. The real part of 
the discrete eigenvalues is reported versus the coupling strength g. 


Let us now compare with the stability of the asynchronous 
solution in the Kuramoto model. Below criticality, in both 
cases the probability distribution is marginally stable (see 
[36] for the first such analysis in the Kuramoto setup), the ma¬ 
jor difference being the presence of a discrete stable spectral 
component in our model. Noteworthy, in spite of the marginal 
stability of the probability density, the order parameter (see 
the next section for its definition) relaxes exponentially in the 
Kuramoto model. This is a manifestation of the so-called 
Landau damping [37]. Only recently this “inconsistency” has 
been fully resolved, by understanding that different classes of 
functions may be considered in the stability analysis [38^0]. 
We do not know how much of such studies carry over to the 
present setup: this is an open problem. 

At criticality, a pair of complex eigenvalues with a posi¬ 
tive real part is born: this is at variance with the standard Ku¬ 
ramoto model, where the newly appearing eigenvalue is real. 
There is, instead, an analogy with the Kuramoto model with 
delay [41, 42], where periodic oscillations arise. Here, how¬ 
ever, above threshold, the probability density rather than os¬ 
cillating periodically, behaves irregularly, as discussed in the 
following sections. 


A. Order parameters 

In order to study the transition, it is necessary to identify a 
suitable order parameter. One cannot directly use the phase 
(p to define the Kuramoto order parameter [25], since in the 
asynchronous regime, such phases do not advance homoge¬ 
neously in time, i.e. they are not proper phases. This can be 
accomplished by introducing the new variable d 

dd ^ t{<P,uj) 

dcP T{l,uj,Eo) ^ ’ 

that is basically equivalent to the elapsed time (apart from 
a scaling factor) and thus advances uniformly by definition. 
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With reference to this new phase, the local dynamics is de¬ 
scribed by the equation 


z? = w - gr{'d){E{t) - Eo) , 
where f (t?) is the effective PRC 


m = 


UJ - gEoT{(/){'&)) 


( 12 ) 


(13) 


'0{4>) is obtained by solving Eq. (11), and a) = l/r(l, oj, Eg), 
is the effective frequency. As it is understood from its defi¬ 
nition depends both on g and the bare frequency. The 
dependence of on ^ is reported in Fig. 3 for the maximal 
and minimal frequencies at g = 0.8. The particular transfor¬ 
mation from the phase </> to the new effective phase i? is shown 
in the appendix B. 



FIG. 3. for g = 0.8 and two different frequencies: uj = 0.8 
(dashed curve) and a; = 2 (solid curve). 


Once a proper phase has been identified, a meaningful 
Kuramoto order parameter can be defined. 






Since i? = 0 in the asynchronous regime, it can be safely used 
to identify the onset of even weak forms of synchronization. 

Given the large amount of transformations needed to de¬ 
termine R, and because of the relationship with neural net¬ 
works, we have often considered a second order parameter, a 
smoothened version Y{t) of the field E[t). In a finite sys¬ 
tem, E{t) is just a collection of (5-pulses. It is therefore more 
convenient to investigate 

Y = -^Y -F E[t) . 


We have selected 7 = 5. In the asynchronous regime, the ac¬ 
tivity is constant, i.e. Yq = Eoj^. Above the transition that 
we are going to discuss, the activity starts oscillating in time, 
so that it is convenient to introduce the temporal standard de¬ 
viation 


= V{y^) - {y? , 


where the angular brackets denote a time average. It will be 
also useful to look at the fluctuation an of the Kuramoto order 
parameter R, as this indicator allows identifying the regimes 
where the degree of synchronizarion oscillates in time. 

We end this theoretical section by briefly commenting on 
the weak-coupling, low-disorder limit. The effective PRC T 
depends on the frequency of the oscillator (and on the cou¬ 
pling strength). The resulting curves for uJmin and uJmax (and 
g = 0.8) are reported in Fig. 1: see the dashed lines. In 
the small-disorder limit, one can neglect such a dependence. 
By then following. Ref. [3], we expect the model to become 
equivalent to the Kuramoto-Daido model 
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The above equation makes one difference with the Kuramoto 
model transparent: the sinusoidal coupling function is re¬ 
placed by the more structured function F (see Fig. 1). This is 
not, however, the major source of differences, since the simu¬ 
lations show that the weak-disorder assumption is not appro¬ 
priate to reproduce the scenario discussed in this paper. 


IV. MACROSCOPIC DYNAMICS 


The most appropriate control parameter to study the onset 
of collective dynamics is the coupling strength g. In Fig. 4a,b, 
it is used to parametrize the dependence of the Kuramoto or¬ 
der parameter R, its temporal standard deviation an, and the 
standard deviation of the activity field Y. 

Each data point is based on a simulation over 500 time units 
after a transient of 50 time units. All, but the red curve, have 
been obtained by increasing the coupling strength g stepwise, 
using the final condition for a given g value as the initial con¬ 
dition for the next one. The statical uncertainty (represented 
by the error bars) has been estimated by dividing the standard 
deviation of the Kuramoto order parameter by the square root 
of the number Ne of effectively independent time intervals 
and Ne has been in turn determined as the ratio between the 
total length of the time series by the decay time of the auto¬ 
correlation function. A long correlation time for p = 1 (se 250 
units) causes the large error compared to smaller {g = 0 . 8 ) 
and larger coupling {g = 1.3). 

The simulations performed for three different system sizes 
(4000, 16000 and 64000 units) reveal the existence of a crit¬ 
ical p-value above which R grows from zero, as in the usual 
Kuramoto setup and that this value is in good agreement with 
gc as estimated from the linear stability analysis discussed in 
the previous section. At variance with the Kuramoto model, 
here, the standard deviation an is larger than zero, meaning 
that the degree of synchronization oscillates in time already 
slightly above threshold. Finally, ay exhibits the same behav¬ 
ior as an, confirming that the transition is accompanied by the 
onset of macroscopic oscillations. 

Finally, the red curve tracks the mean-field obtained by de¬ 
creasing g. The difference observed in the critical region with 
respect to the previous curves (obtained by increasing g) sug- 
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FIG. 4. Phase diagram; Dependence of the Kuramoto order parame¬ 
ter R (a), its standard deviation aR, the standard deviation ay of the 
activity field Y (b), and the ten largest Lyapunov exponents versus 
the coupling strength g. The curves in a and b have been obtained 
for N = 4000 (black), N = 16000 (green) and N = 64000 (blue), 
upon increasing the coupling g. The red curve (a and b) is a con¬ 
tinuation by decreasing g with N — 64000. or and ay are repre¬ 
sented with dashed and dotted lines, respectively (b). The vertical 
lines mark the critical point g^ ~ 0.72, where the linear stability of 
the asynchronous state is lost (see section III). The ten largest global 
Lyapunov exponents in panel c are almost indistinguishable. The 
simulations have been performed for N = 64000, increasing g. 


gests the possible co-existence of an asynchronous with a par¬ 
tially synchronised regime. Since, however, no jump is ob¬ 
served in the simulations performed by increasing g, it is rea¬ 
sonable to conclude that the bifurcation is “supercritical” and 
thus to attribute such deviations to the finite sweeping time. 
Anyway, since the main goal of this work is to characterise the 
behavior above threshold, we have preferred to focus our ef¬ 
forts on larger p-values, where the asymptotic regime is much 
less dependent on the selection of the initial condition. 

A first qualitative instance of the collective dynamics can 
be appreciated in Fig. 5, where R{t) and Y(t) are plotted for 
g = 1.3 showing that the evolution is more complex than just 
periodic. 

A more quantitative characterization of the collective tem¬ 
poral behavior can be obtained by looking at power spectra. 
The square amplitude of the Fourier transform of Y (t) is re¬ 
ported in Fig. 6 for two different coupling strengths. 

There we see that the spectra possess quite a rich structure, 
being neither trivially broad-band, nor just revealing a peri¬ 
odic behavior (especially for g = 0.8). A closer look at the 
width of the various peaks upon increasing the network size 
reveals that they do not decrease. Simulations performed for 
different realization of the bare frequencies (data not shown) 
indicate that the results are almost independent, especially for 
the larger system sizes. Altogether, we are thus led to con¬ 
clude that the stochastic-like dynamics is not due to finite-size 
effects, but intrinsic of the thermodynamic limit. 


0 5 10 15 



FIG. 5. Time evaluation of Y (a) and R (b) for g — 1.3 and system 
sizes: N = 4000. 
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FIG. 6. Power spectra of T for p = 0.8 (a), and 1.3 (b) in each case 
for: N = 4000 (black), N = 16000 (red) and N = 64000 (green). 
The spectra are obtained by transforming time series of 819.175 time 
units, sampled every 0.025 units and averaged over 50 different real¬ 
izations. 


According to the theory of nonlinear dynamical systems, 
it is well known that an irregular evolution may well be the 
manifestation of low-dimensional deterministic chaos. Can it 
be the case here? In order to clarify the point, it is natural 
to investigate the behavior of the activity field by perform¬ 
ing a nonlinear time series analysis, to determine its fractal 
dimension. Given a time series Y{tn), sampled at equally 
spaced times (At = — tn = 0.025), one starts embed¬ 

ding the series into a space of dimension m, by building vec¬ 
tors of the type [Y (tn),Y (t^+s), ■ ■ ■ ,Y (f„+(m_i)s)], where 
s is suitably selected. As often done, we have chosen s, so 
that sAt is close to the first minimum of the autocorrelation 
of Y(t) (s = 5, in our case). 

The fractal dimension has then been estimated by using the 
nearest-neighbour method [43], as it suffers of less fluctua- 
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tions in the region of small distances. Given a generic time 
series, reference points are randomly selected {N^. = 10® 
in our case). Each of them is compared with an increasing 
number n of randomly selected measurement points (the other 
points in the time series - up to a maximum Nm = 16 • 10®), 
monitoring the distance em{k,n) of the fc-th neighbour (the 
distance is herein estimated using the maximum norm), for 
different values of the embedding dimension m and k. A well 
established theory [43], implies that for large n 

{\nem{k,n)) ^ , 


where the angular brackets denote the average over the refer¬ 
ence points, while Dg is the (effective) information dimen¬ 
sion. In order to make the dependence of De on the res¬ 
olution Em transparent, we have modified the standard ap¬ 
proach. Once interpreted the logarithmic derivative of n as 
a resolution-dependent dimension. 


(^m) 


dlnn 

d{liiem) 


we have plotted it versus (Em) itself, interpreted as an inde¬ 
pendent variable. In fact, we have verified that DeiEm) takes 
the same value, irrespective of the way Em has been deter¬ 
mined (i.e. independently of the k value). The only differ¬ 
ences are that larger fc-values yield smaller statistical fluctu¬ 
ations, but are confined to larger distances. A good compro¬ 
mise has been obtained by gluing together the data obtained 
for the largest k value (30) with the data obtained for the 
smallest distances and a lower-order neighbour (4th one). 



not even infinite dimensional). Additionally, one can also ap¬ 
preciate a small shift to smaller scales of the curves obtained 
for N = 64000. In itself, this is the indication of finite-size 
effects. If the shift continues as such by further increasing 
m, it would mean that part of the high-dimensionality is just 
a consequence of statistical fluctuations which disappear in 
the thermodynamic limit. We are more inclined to attribute 
such discrepancies to another type of finite-size effect: a non 
perfect equivalence among the various realizations of the fre¬ 
quency distributions. We have, in fact, observed that different 
clusters may temporarily form during the evolution, especially 
in the interval g G [gc, 1-2] (see a more detailed discussion in 
section V). 

For p = 1.3 the convergence to the thermodynamic limit is 
more clear. In Fig. 7b the agreement among the different net¬ 
work sizes is compelling over a wide range, suggesting thus 
that the statistical fluctuations do not affect the dimension es¬ 
timates. We have, also, double-checked the results by com¬ 
puting the correlation dimension with the TISEAN package 
[44]: a rather similar pattern emerges (data not shown). 

Finally, we have investigated the degree of (in)stability of 
the dynamics, computing the first 10 Fyapunov exponents Aj. 
We have followed the approach described in [45], which con¬ 
sists in formally interpreting the time evolution as a series 
of discrete-time maps from one to the next spike emission. 
The results are plotted in Fig. 4c upon varying the coupling 
strength. There we see that the dynamics is always stable (no¬ 
tice that the zero Fyapunov exponent, always present in a non¬ 
constant autonomous dynamics is automatically discarded). 
For a vanishing g, all the Fyapunov exponents converge to 
zero: this is obvious, since in this limit all the oscillators are 
uncoupled. Much less trivial is that the Fyapunov exponents 
are all negatives in spite of a dynamics that may even be col¬ 
lectively irregular. This manifestation of stable chaos strongly 
suggests that the connection between different levels of de¬ 
scriptions (micro vs. macroscopic) of a given model is weak, 
if any. 


V. MICROSCOPIC DYNAMICS 


FIG. 7. Effective dimension De as a function of the resolution Em for 
g — 0.8 and g = 1.3 in panel a and b, respectively. The system size 
is A = 4000 (black), N = 16000 (red) and N — 64000 (green). 
The different symbols belong to different embedding dimensions m 
marked in the figure. The curves for the same embedding dimension 
group together with a similar slope irrespective of the system size. 


The results for g = 0.8 are reported in Fig. 7a and differ¬ 
ent ensemble sizes {N = 4000, 16000, and 64000). The first 
point to notice is that for the lowest embedding dimension 
(m = 4), De nicely converges to 4 upon decreasing e. This 
clearly implies that the dimension of the collective motion is 
at least 4, i.e. one needs at least four variables to character¬ 
ize such a behavior. Furthermore, the curves obtained for the 
larger m values, reveal an increase, each possibly hinting to 
m, thus suggesting that the dynamics is high-dimensional (if 


In this section we try to shed light on the collective dynam¬ 
ics by analysing the behavior of the single neurons. We start 
noticing that the coupling modifies the firing rate of the neu¬ 
rons. This can be appreciated in Fig. 8a, where the effective 
(average) frequency Cj is reported for the coupling strength 
p = 1. The dashed line corresponds to the bare frequency of 
each neuron. Almost everywhere Cj is smaller than the bare 
frequency w. This is a consequence of the fact that, although 
the PRC was chosen to be symmetric around zero, this is no 
longer true for the effective PRC (see Fig. 1). The most inter¬ 
esting feature to notice is, however, the staircase structure of 
Cj{ijj) with flat plateaus which correspond to clusters of mutu¬ 
ally synchronised neurons: the synchronization does not mean 
a perfect phase locking but that the phase differences never 
become larger than 27r. 

One way to characterize the irregularity of the single¬ 
neuron activity is through its coefficient of variation (CV), i.e. 
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the standard deviation of the the inter-spike interval rescaled 
to its average value. In Fig. 8b, we see that the CV allows 
identifying synchronized clusters as those frequency intervals 
where the fluctuations are significantly smaller. Furthermore, 
distinct lines can be recognized inside some clusters: they cor¬ 
respond to different locked states [46] and are a manifestation 
of the multistability that is in fact seen also at the macroscopic 
level. On a more quantitative level, the neural dynamics is not 
significantly irregular if compared, for instance, to the true 
brain activity in the resting state. It should be, however, kept 
in mind that in our toy model, the only source of disorder is 
the distribution of bare frequencies; no disorder has been as¬ 
sumed in the synaptic connections. 

Additional information can be extracted by assuming that 
the self-determined activity field E{t) is externally given, so 
that each neuron can be interpreted as a forced dynamical sys¬ 
tem. In this way, it is natural to compute the (conditional) 
Lyapunov exponent Ai. In Fig. Ic, one can observe a sce¬ 
nario that is qualitatively different from what observed in the 
Kuramoto model: all A^’s are negative, including those of the 
neurons outside the flat plateaus. We come back to this point 
ahead in this section. 

A partially different scenario is found for p = 1.3 (see 
Fig. 8d-f). First of all any sign of multistability has disap¬ 
peared and all plateaus as well. The initial high peak of the 
CV (panel e) is due to the fact that now the neurons with 
the lowest bare frequency do not spike at all, having under¬ 
gone a kind of oscillation-death. Their CV is obviously equal 
to zero. As a result, the first erraticaly spiking neurons are 
characterized by long interspike intervals and are obviously 
accompanied by large fluctuations. The CV of these rarely br¬ 
ing neurons is as large as 1.75 (not shown in Fig. 8e with the 
present scale). The kind of oscillation-death phenomena could 
be seen as an inhomogeneous limit cycle (IHLC) because non 
spiking neurons are not trapped at a steady state, rather mov¬ 
ing back and forth according to their bare frequency and the 
global pulses, never reaching the threshold [47]. Hence we 
observe two groups of neurons, the quiet neurons but still be¬ 
ing sub-thresholdly active and the bring neurons. The situa¬ 
tion is a more complicated than usual IHLC, because of the 
many (in the thermodynamic limit infinite many) frequencies 
reflected in the power spectra (Fig. 6). Moreover, the dynam¬ 
ics of each oscillator is rather stable. 

Altogether, the microscopic analysis confirms that the mi¬ 
croscopic behavior is linearly stable: each neuron is synchro¬ 
nized with the self-generated mean-field E{t) and yet an ir¬ 
regular dynamics is self-sustained. This can be understood 
by noticing that in many frequency ranges, the effective fre¬ 
quency is a strictly monotonous function of tu. This means 
that, even though each neuron synchronizes with the held E, 
the parameter (bare frequency) mismatch induces a qualita¬ 
tively different response: such qualitative differences are then 
responsible for the maintenance of the self-generated irregu¬ 
larity. In a more technical way: the response of a phase os¬ 
cillator is not structurally stable: one can slightly modify its 
frequency and still observe significant changes (phase-slips). 
A more physical (although still qualitative) way to under¬ 
stand the phenomenon is as follows: one can see each phase- 




FIG. 8. Effective frequency Cj (panels a and d), CV of ui (panels b 
and e) and conditional Lyapunov exponents Ai (panels c and f) for 
the different oscillators, for g = 1 (panels a-c) and g = 1.3 (panels 
d-f), all for N — 4000. The red dashed line in panels a and d show 
the bare frequency as a reference for the effective frequency shown 
with solid black lines. 


oscillator as a particle moving in a potential with an inclina¬ 
tion that depends on its bare frequency and the mean-field E. 
When two particles with slightly different w are followed, it 
may happen that one of them is blocked in a (shallow) mini¬ 
mum which is absent for the other. 

One can learn a bit more about the dependence on uj by 
comparing pairs of consecutive oscillators (consecutive in the 
space of bare frequencies). This can be done, by monitor¬ 
ing phase slips, i.e. the time instants when the phase differ¬ 
ence becomes larger (smaller) than 1/2 (—1/2). Since it is 
possible that the phase difference may oscillate around 1/2 
(—1/2) yielding long sequences of positive and negative slips, 
we have chosen to record only those events where two or more 
consecutive positive (negative) flips are observed. 

In Fig. 9a we see that for 5 = 1 there exist totally empty 
bands: they correspond to the previously mentioned synchro¬ 
nization areas. In panel b (which corresponds to g = 1.3) 
both forward and backward slips are simultaneously present. 
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FIG. 9. Phase slips for g = 1 (panel a) and g = 1.3 (panel b) for 
N = 4000. Black circles correspond to forward slips; red crosses to 
backward slips, occurring only in panel b. 


One can see that the phase slips happen on a much longer 
time scale than the interspike intervals. Further, the bands 
with sparse phase slips observed for g = 1.3 are reminiscent 
of the synchronization bands found for g = 1.0: there, phase 
slip events are rare and erratic. The totally empty band at the 
bottom frequencies for g = 1.3 corresponds to the non-hring 
neurons. 

In the thermodynamic limit, the most appropriate way to 
characterize the collective dynamics is by monitoring the 
probability distribution Q{<j),Lo,t) introduced in Sec. III. In 
Fig. 10 we give an idea of the way it looks like below and 
above threshold at some randomly chosen time. In panel a) 
one can recognize a reasonably smooth distribution. In fact, 
for g = 0.5, we are in the asynchronous regime and thereby 
expect a smooth distribution of the phases <j) [48] Such a dis¬ 
tribution looses stability above gc- 



FIG. 10. Snapshots of the probability density t) for g = 0.5 

(panel a) and <; = 1.3 (panel b) and N = 4000. 


In fact, for g = 1.3 we see a rather different structure (see 
panel b) characterized by an alternation of highly dense and 
widely spread regions. It is clear that even the plain integra¬ 
tion of the equations (4,5) is a highly nontrivial task, not to 
speak of the development of approximate analytical schemes. 


VI. DISCUSSION AND OPEN PROBLEMS 


In this paper we have analysed an ensemble of pulse- 
coupled oscillators characterized by a distribution of bare fre¬ 
quencies and coupled through a homogeneous mean held. Al¬ 
though the setup is reminiscent of the Kuramoto model, the 
collective dynamics is much richer and accompanied by a lin¬ 
early stable microscopic dynamics. 

A linear stability analysis of the asynchronous regime al¬ 
lows identifying the transition point beyond which a com¬ 
plex form of synchronization sets in. A numerical analysis 
of a properly dehned Kuramoto order parameter R and of 
the smoothed activity held Y reveals that they not only huc- 
tuate in time, but their behavior involves a large (possibly 
inhnite) number of degrees of freedom. This indicates that 
even in “simple” mean-held models, such as the one investi¬ 
gated in this paper, the coarse-grained activity of an ensem¬ 
ble of phase-oscillators cannot be reduced to the evolution 
of one or a few variables, such as the bring rate and related 
observables. In principle, nothing prevents a population of 
phase-oscillators to self-sustain a macroscopic irregular dy¬ 
namics: the corresponding evolution equation is indeed a non¬ 
linear functional equation (see Eqs. (4,5)), which operates in 
an inhnite-dimensional phase-space. It is, however, unclear 
under which conditions many degrees of freedom can be si¬ 
multaneously active. In order to make further progress, it will 
be necessary to hnd suitable approximations of the probabil¬ 
ity density Q{<j),Lo,t): this task seems to require clever ideas 
on the way to expand One question is particu¬ 

larly relevant: whether the dynamics is born high dimensional 
from the very beginning (such as in models of balanced states 
[49, 50]) or the complexity increases by undergoing a series of 
consecutive bifurcations. The numerical analysis in the vicin¬ 
ity of the critical point is affected by too strong hnite-size cor¬ 
rections to be able to draw any conclusion. 

Another open point is the generality of this scenario. Sev¬ 
eral preliminary simulations performed with various choices 
of the PRC reveal that it is quite robust, although the presence 
of a relatively steep branch seems to be a necessary condition. 
This is not too serious a limitation, as it naturally appears in 
systems characterized by a slow-fast dynamics (see the dis¬ 
cussion in [51]). It might be, however, worth to assume a dif¬ 
ferent PRC shape to enable deeper analytical studies. We have 
indeed derived a very general equation for the loss of stability 
of the asynchronous state; if one could go beyond, including 
the most relevant nonlinear terms, it should be possible to de¬ 
cide how many degrees of freedom are switched on. 

Our numerical studies suggest that the transition disappears 
when the distribution of frequencies is narrow enough, but this 
is by no means a proof: understanding whether it is strictly 
necessary to go beyond the weak-coupling, weak-disorder 
limit is another point that will be worth exploring. 

Another intriguing property of the collective dynamics dis¬ 
cussed in this paper is the presence of a spectrum of negative 
Lyapunov exponents. This means that it is a manifestation of 
stable chaos [29]. Within the context of computational neuro¬ 
science, the stability of the microscopic trajectories suggests 
that this model is a good candidate for performing computa- 
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tional tasks. It will be worth to explore this opportunity by 
studying the response of this type of networks to different 
classes of external stimuli. 
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Appendix A: Stability analysis 

The equation (9) for u can be rewritten as 

du _ {9^\4 ')Eo-^J) u +g[T'{(j))Qo + T{(j))^]z 
d 4 > uj-gT{(j))Eo 


which has the structure 
du 

— = (A - pt)m + gP{uj)C{(j), ui)z 

where 

gT'{4’)Eo 

^-9^{(t>)Eo 


(Al) 


(A2) 


C'( 0 ,cu) = 


ccr'(^) 


(A3) 


T{uj,Eo){LU-gT{^)Eo)^ 

while T is defined in (7). The general solution of Eq. (Al) is 

u{^, uj) = [u(o, uj)+ (A4) 

Jo 


where 


ri’ 

E{'ip,uj)= / d77A(p)=log 
Jo 


^ - 9^{0)Eo 

^-9^{'J’)Eo 


(A5) 


notice that F{1) = F{Q) = 0, while T{ip,uj) is defined in 
Eq. (7). One can therefore write the solution u{(j), w) as 

UJ - gr[4>)Eo 


zgP{uj) 




UJ - gT{0)Eo 


where 


V^{(j),uj) = 


r(i 


(l>w) Jo 


d'ip- 


r'(V')e 




- 9^{J’)Eof 


(A6) 


(A7) 


By imposing the periodicity condition u(l,a;) = u{0,uj) one 
obtains 


=9Z 


P{uj)Vf,{l,uj){l,uj) 
(eMT(i,<.) _i)(^_gr( 0 )E;o) 


(A8) 


By finally, recalling the definition (10) z, it is possible to 
rewrite Eq. (A8) as the eigenvalue equation 


l + 5 r(l) J Qo{1,uj)(Ej = g J dujJ^^^^^J^ (A9) 

Appendix B: Transformation to the appropriate Kuramoto like 
phase 9 


We solve the original system Eq. (1) with the PRC Eq. (2) 
in terms of the non-homogeneously advancing phase (jj. At 
each time point we calculate the Kuramoto order parameter 
R we convert the phase (jj into the Kuramoto like phase d ac¬ 
cording to Eqs. (12) and (13). Eor simplicity and readability 
we introduce new constants similar those defining the piece- 
wise linear PRC Eq. (2): *8oi = -^EqBqi, $02 = -^EqBq 2 , 
®03 = jfEoBQs, hi = j^Eobi, and b 2 = jjEah 2 . Hence the 
Kuramoto like phase is 


t?((^(f)) = < 


£L In 

bi w —'Boi — 

,9 _^ — 2?02 + b 2 <i>i 

^ b2 cj —t8o2-t-b2(^(t) 

-1- J2. In ‘^-‘So3-bl0T- 

Fr -I- bi u,-<&03-bi<l>(t) 


if 0 < Ip < 4'1 
if Pi <P<Pr 

if pr < P < I , 


with the field Eq as stated in Eq. (8) and the effective fre¬ 
quency defined as the inverse interspike interval, i.e. uj = 
1/T{l,uj, Eq) (Eq. (7)). The interspike interval can be ex¬ 
pressed explicitly for the give PRC: 


T{1,uj,Eq) = — In- - -— 

bi w-Soi-bn; 


b2pi 


1 w — *Bo 2 

b2 UJ - ®02 + i>2pr 

1 . UJ — ^03 — bipr 
bi w - $03 - fat 


As shown in Pig. 1, the transitions 1)1 and dr in the effective 
PRC r depend on the bare frequency uj according to: 


uj UJ - $01 

— In- 

bi uj-^Qi-bipi 

_ UJ - $02 -I- b 2 ^i 

^ b 2 ^ w - $02 + b2pr 
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